Plot standard errors (SE) from Local h2 estimation of Cross-Tissue and Tissue-Wide expression in GTex

library(ggplot2)
library(reshape2)
library(dplyr)
library(tidyr)
library(GGally)
"%&%" = function(a,b) paste(a,b,sep="")
#my.dir <- '/Users/heather/Dropbox/cross-tissue/expArch_DGN-FHS_results/GTEx-cross-tissue/'
my.dir <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/'
tislist <- scan(my.dir %&% 'GTEx_nine_tissues','c')

ct <- read.table(my.dir %&% 'GTEx.ranef.cross-tissue.h2_marginal.local_2015-03-23.txt',header=T)
nlist<-ct[1,2]

allse <- ct %>% select(ensid,se)
colnames(allse) <- c('ensid','CrossTissue')

for(i in 1:length(tislist)){
  tis <- tislist[i]
  tisdata <- read.table(my.dir %&% 'GTEx.tissue-wide.h2_' %&% tis %&% '_marginal.local_2015-03-24.txt',header=T)
  nlist <- c(nlist,tisdata[1,2])
  tisse <- tisdata %>% select(ensid,se)
  allse <- inner_join(allse,tisse,by='ensid')
  tis <- gsub('(','',tis,fixed=T) #replace any parentheses in tis string
  tis <- gsub(')','',tis,fixed=T)
  tis <- gsub('-','',tis)
  colnames(allse)[2+i] <- tis
}

alllist<-c("Cross-Tissue",tislist)
samplesize<-data.frame(cbind(alllist,nlist))
colnames(samplesize)<-c("Tissue","n")
arrange(samplesize,desc(n))
##                       Tissue   n
## 1               Cross-Tissue 450
## 2            Muscle-Skeletal 361
## 3                 WholeBlood 339
## 4  Skin-SunExposed(Lowerleg) 303
## 5       Adipose-Subcutaneous 298
## 6              Artery-Tibial 285
## 7                       Lung 279
## 8                    Thyroid 279
## 9               Nerve-Tibial 256
## 10       Heart-LeftVentricle 190
dim(allse)
## [1] 17007    11
a<-gather(allse,"CrossTissue","Tissue",3:11)
colnames(a) <- c('ensid','CrossTissue','TissueName','Tissue')
ggplot(a,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(cex=0.7) + geom_abline(intercept=0, slope=1,color='red') + geom_smooth(method = "lm") + ylab('Cross-Tissue SE') + xlab('Tissue-Wide SE') +  ggtitle("GTEx Compare SE from Local h2 estimation")

Plot Local h2 estimates of Cross-Tissue and Tissue-Wide expression in GTex

allh2 <- ct %>% select(ensid,h2)
colnames(allh2) <- c('ensid','CrossTissue')

for(i in 1:length(tislist)){
  tis <- tislist[i]
  tisdata <- read.table(my.dir %&% 'GTEx.tissue-wide.h2_' %&% tis %&% '_marginal.local_2015-03-24.txt',header=T)
  tish2 <- tisdata %>% select(ensid,h2)
  allh2 <- inner_join(allh2,tish2,by='ensid')
  tis <- gsub('(','',tis,fixed=T) #replace any parentheses in tis string
  tis <- gsub(')','',tis,fixed=T)
  tis <- gsub('-','',tis)
  colnames(allh2)[2+i] <- tis
}

dim(allh2)
## [1] 17007    11
b<-gather(allh2,"CrossTissue","Tissue",3:11)
colnames(b) <- c('ensid','CrossTissue','TissueName','Tissue')
ggplot(b,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(cex=0.7) + geom_abline(intercept=0, slope=1,color='red') + geom_smooth(method = "lm") + ylab('Cross-Tissue h2') + xlab('Tissue-Wide h2') +  ggtitle("GTEx Local h2 estimates")

Plot histograms of h2 and se

histh2<-gather(allh2,"ensid","Tissue",2:11)
colnames(histh2) <- c('ensid','Tissue','h2')
ggplot(histh2, aes(x = h2, fill = Tissue)) + geom_histogram() + facet_wrap(~Tissue,ncol=2) +  ggtitle("GTEx Local h2") + guides(fill=FALSE)

ggplot(histh2, aes(x = h2, fill = Tissue)) + geom_histogram() + facet_wrap(~Tissue,ncol=2) + coord_cartesian(xlim=c(0.2,1),ylim=c(0,300)) + ggtitle("GTEx Local h2 > 0.20") + guides(fill=FALSE)

ggplot(histh2, aes(x = h2, fill = Tissue, color = Tissue)) + geom_density(alpha=0.3) + ggtitle("GTEx Local h2")

ggplot(histh2, aes(x = h2, fill = Tissue, color = Tissue)) + geom_density(alpha=0.3) + coord_cartesian(ylim=c(0,15)) + ggtitle("GTEx Local h2 ZOOM")

histse<-gather(allse,"ensid","Tissue",2:11)
colnames(histse) <- c('ensid','Tissue','se')
ggplot(histse, aes(x = se, fill = Tissue, color = Tissue)) + geom_density(alpha=0.3)

Do same with tissue-specific h2 and se

Plot standard errors (SE) from Local h2 estimation of Cross-Tissue and Tissue-Specific expression in GTex

my.dir <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/'
tislist <- scan(my.dir %&% 'GTEx_nine_tissues_spaces','c',sep="\n")

ct <- read.table(my.dir %&% 'GTEx.ranef.cross-tissue.h2_marginal.local_2015-03-23.txt',header=T)
nlist<-ct[1,2]

ts_allse <- ct %>% select(ensid,se)
colnames(ts_allse) <- c('ensid','CrossTissue')

for(i in 1:length(tislist)){
  tis <- tislist[i]
  tisdata <- read.table(my.dir %&% 'GTEx.resid.tissue-specific.h2_' %&% tis %&% '_marginal.local_2015-03-23.txt',header=T,sep="\t")
  nlist <- c(nlist,tisdata[1,2])
  tisse <- tisdata %>% select(ensid,se)
  ts_allse <- inner_join(ts_allse,tisse,by='ensid')
  tis <- gsub(' ','',tis) #replace spaces in tis string
  tis <- gsub('(','',tis,fixed=T) #replace any parentheses in tis string
  tis <- gsub(')','',tis,fixed=T)
  tis <- gsub('-','',tis)
  colnames(ts_allse)[2+i] <- tis
}

alllist<-c("Cross-Tissue",tislist)
samplesize<-data.frame(cbind(alllist,nlist))
colnames(samplesize)<-c("Tissue","n")
arrange(samplesize,desc(n))
##                            Tissue   n
## 1                    Cross-Tissue 450
## 2               Muscle - Skeletal 361
## 3                     Whole Blood 339
## 4  Skin - Sun Exposed (Lower leg) 303
## 5          Adipose - Subcutaneous 298
## 6                 Artery - Tibial 285
## 7                            Lung 279
## 8                         Thyroid 279
## 9                  Nerve - Tibial 256
## 10         Heart - Left Ventricle 190
dim(ts_allse)
## [1] 17022    11
ts_a<-gather(ts_allse,"CrossTissue","Tissue",3:11)
colnames(ts_a) <- c('ensid','CrossTissue','TissueName','Tissue')
ggplot(ts_a,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(cex=0.7) + geom_abline(intercept=0, slope=1,color='red') + geom_smooth(method = "lm") + ylab('Cross-Tissue SE') + xlab('Tissue-Specific SE') +  ggtitle("GTEx Compare SE from Local h2 estimation")

Plot Local h2 estimates of Cross-Tissue and Tissue-Specific expression in GTex

ts_allh2 <- ct %>% select(ensid,h2)
colnames(ts_allh2) <- c('ensid','CrossTissue')

for(i in 1:length(tislist)){
  tis <- tislist[i]
  tisdata <- read.table(my.dir %&% 'GTEx.resid.tissue-specific.h2_' %&% tis %&% '_marginal.local_2015-03-23.txt',header=T,sep="\t")
  tish2 <- tisdata %>% select(ensid,h2)
  ts_allh2 <- inner_join(ts_allh2,tish2,by='ensid')
  tis <- gsub(' ','',tis) #replace spaces in tis string
  tis <- gsub('(','',tis,fixed=T) #replace any parentheses in tis string
  tis <- gsub(')','',tis,fixed=T)
  tis <- gsub('-','',tis)
  colnames(ts_allh2)[2+i] <- tis
}

dim(ts_allh2)
## [1] 17022    11
ts_b<-gather(ts_allh2,"CrossTissue","Tissue",3:11)
colnames(ts_b) <- c('ensid','CrossTissue','TissueName','Tissue')
ggplot(ts_b,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(cex=0.7) + geom_abline(intercept=0, slope=1,color='red') + geom_smooth(method = "lm") + ylab('Cross-Tissue h2') + xlab('Tissue-Specific h2') +  ggtitle("GTEx Local h2 estimates")

Plot histograms of h2 and se Cross-Tissue and Tissue-Specific

histh2<-gather(ts_allh2,"ensid","Tissue",2:11)
colnames(histh2) <- c('ensid','Tissue','h2')
ggplot(histh2, aes(x = h2, fill = Tissue)) + geom_histogram() + facet_wrap(~Tissue,ncol=2) +  ggtitle("GTEx Local h2") + guides(fill=FALSE)

ggplot(histh2, aes(x = h2, fill = Tissue)) + geom_histogram() + facet_wrap(~Tissue,ncol=2) + coord_cartesian(xlim=c(0.2,1),ylim=c(0,300)) + ggtitle("GTEx Local h2 > 0.20") + guides(fill=FALSE)

ggplot(histh2, aes(x = h2, fill = Tissue, color = Tissue)) + geom_density(alpha=0.3) + ggtitle("GTEx Local h2")

ggplot(histh2, aes(x = h2, fill = Tissue, color = Tissue)) + geom_density(alpha=0.3) + coord_cartesian(ylim=c(0,15)) + ggtitle("GTEx Local h2 ZOOM")

histse<-gather(ts_allse,"ensid","Tissue",2:11)
colnames(histse) <- c('ensid','Tissue','se')
ggplot(histse, aes(x = se, fill = Tissue, color = Tissue)) + geom_density(alpha=0.3)

se3<-inner_join(allse,ts_allse,by="ensid")
h23<-inner_join(allh2,ts_allh2,by="ensid")

seTW <- select(se3,1:11)
colnames(seTW)<-gsub('\\.x','',colnames(seTW))
gseTW <- gather(seTW,"Tissue","TissueWide",3:11)
dim(gseTW)
## [1] 153063      4
seTS <- select(se3,ensid,12:21)
colnames(seTS)<-gsub('\\.y','',colnames(seTS))
gseTS <- gather(seTS,"TissueName","TissueSpecific",3:11)
dim(gseTS)
## [1] 153063      4
gseALL <- cbind(gseTW,gseTS)


h2TW <- select(h23,1:11)
colnames(h2TW)<-gsub('\\.x','',colnames(h2TW))
gh2TW <- gather(h2TW,"Tissue","TissueWide",3:11)
dim(gh2TW)
## [1] 153063      4
h2TS <- select(h23,ensid,12:21)
colnames(h2TS)<-gsub('\\.y','',colnames(h2TS))
gh2TS <- gather(h2TS,"TissueName","TissueSpecific",3:11)
dim(gh2TS)
## [1] 153063      4
gh2ALL <- cbind(gh2TW,gh2TS)

ggplot(gh2ALL,aes(x=TissueSpecific,y=TissueWide)) +facet_wrap(~Tissue,scales="fixed",ncol=3) + geom_point(cex=0.7) + geom_abline(intercept=0, slope=1,color='red') + geom_smooth(method = "lm") + ylab('Tissue-Wide h2') + xlab('Tissue-Specific h2') +  ggtitle("GTEx Local h2 estimates")

ggplot(gseALL,aes(x=TissueSpecific,y=TissueWide)) +facet_wrap(~Tissue,scales="fixed",ncol=3) + geom_point(cex=0.7) + geom_abline(intercept=0, slope=1,color='red') + geom_smooth(method = "lm") + ylab('Tissue-Wide SE') + xlab('Tissue-Specific SE') +  ggtitle("GTEx Compare SE from Local h2 estimation")

ggplot(gh2ALL,aes(x=TissueSpecific+CrossTissue,y=TissueWide)) +facet_wrap(~Tissue,scales="fixed",ncol=3) + geom_point(cex=0.7) + geom_abline(intercept=0, slope=1,color='red') + geom_smooth(method = "lm") + ylab('Tissue-Wide h2') + xlab('Tissue-Specific + Cross-Tissue h2') +  ggtitle("GTEx Local h2 estimates")

###make ordered point+CI h2 plots
gTW_h2 <- gather(h2TW,"Tissue.h2","TissueWide.h2",2:11)
gTW_se <- gather(seTW,"Tissue.se","TissueWide.se",2:11)
gTW_h2_se <- cbind(gTW_h2,gTW_se[,3])
colnames(gTW_h2_se) <- c("ensid","Tissue","h2","se")
ngenes<-nrow(gTW_h2_se)/length(unique(gTW_h2_se$Tissue))
gTW_h2_se <- gTW_h2_se %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) ) %>% arrange(Tissue,h2) %>% mutate(place=rep(1:ngenes,10),nonzeroCI=ymin>0) 
ggplot(gTW_h2_se,aes(x=place,y=h2,ymin=ymin, ymax=ymax, color=nonzeroCI) ) + facet_wrap(~Tissue,ncol=2) + geom_pointrange(col='gray')+geom_point()+ggtitle('Tissue-Wide h2')+ theme(legend.justification=c(0,1), legend.position=c(0,1))+ylab(expression("local h"^2)) + xlab(expression("genes ordered by local h"^2))

gTS_h2 <- gather(h2TS,"Tissue.h2","TissueSpecific.h2",2:11)
gTS_se <- gather(seTS,"Tissue.se","TissueSpecific.se",2:11)
gTS_h2_se <- cbind(gTS_h2,gTS_se[,3])
colnames(gTS_h2_se) <- c("ensid","Tissue","h2","se")
ngenes<-nrow(gTS_h2_se)/length(unique(gTS_h2_se$Tissue))
gTS_h2_se <- gTS_h2_se %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) ) %>% arrange(Tissue,h2) %>% mutate(place=rep(1:ngenes,10),nonzeroCI=ymin>0)
ggplot(gTS_h2_se,aes(x=place,y=h2,ymin=ymin, ymax=ymax, color=nonzeroCI) ) + facet_wrap(~Tissue,ncol=2) + geom_pointrange(col='gray')+geom_point()+ggtitle('Tissue-Specific h2')+ theme(legend.justification=c(0,1), legend.position=c(0,1))+ylab(expression("local h"^2)) + xlab(expression("genes ordered by local h"^2))

##save data for later
write.table(h2TW,file=my.dir %&% "GTEx_Tissue-Wide_local_h2.txt",quote=F,row.names=F)
write.table(seTW,file= my.dir %&% "GTEx_Tissue-Wide_local_se.txt",quote=F,row.names=F)
write.table(h2TS,file=my.dir %&% "GTEx_Tissue-Specific_local_h2.txt",quote=F,row.names=F)
write.table(seTS,file=my.dir %&% "GTEx_Tissue-Specific_local_se.txt",quote=F,row.names=F)

TW<-inner_join(h2TW,seTW,by="ensid")
comboTW <- transmute(TW,EnsemblGeneID=ensid,h2.CrossTissue=CrossTissue.x,se.CrossTissue=CrossTissue.y,h2.AdiposeSubcutaneous=AdiposeSubcutaneous.x,se.AdiposeSubcutaneous=AdiposeSubcutaneous.y,h2.ArteryTibial=ArteryTibial.x,se.ArteryTibial=ArteryTibial.y,h2.HeartLeftVentricle=HeartLeftVentricle.x,se.HeartLeftVentricle=HeartLeftVentricle.y,h2.Lung=Lung.x,se.Lung=Lung.y,h2.MuscleSkeletal=MuscleSkeletal.x,se.MuscleSkeletal=MuscleSkeletal.y,h2.NerveTibial=NerveTibial.x,se.NerveTibial=NerveTibial.y,h2.SkinSunExposedLowerleg=SkinSunExposedLowerleg.x,se.SkinSunExposedLowerleg=SkinSunExposedLowerleg.y,h2.Thyroid=Thyroid.x,se.Thyroid=Thyroid.y,h2.WholeBlood=WholeBlood.x,se.WholeBlood=WholeBlood.y)

TS<-inner_join(h2TS,seTS,by="ensid")
comboTS <- transmute(TS,EnsemblGeneID=ensid,h2.CrossTissue=CrossTissue.x,se.CrossTissue=CrossTissue.y,h2.AdiposeSubcutaneous=AdiposeSubcutaneous.x,se.AdiposeSubcutaneous=AdiposeSubcutaneous.y,h2.ArteryTibial=ArteryTibial.x,se.ArteryTibial=ArteryTibial.y,h2.HeartLeftVentricle=HeartLeftVentricle.x,se.HeartLeftVentricle=HeartLeftVentricle.y,h2.Lung=Lung.x,se.Lung=Lung.y,h2.MuscleSkeletal=MuscleSkeletal.x,se.MuscleSkeletal=MuscleSkeletal.y,h2.NerveTibial=NerveTibial.x,se.NerveTibial=NerveTibial.y,h2.SkinSunExposedLowerleg=SkinSunExposedLowerleg.x,se.SkinSunExposedLowerleg=SkinSunExposedLowerleg.y,h2.Thyroid=Thyroid.x,se.Thyroid=Thyroid.y,h2.WholeBlood=WholeBlood.x,se.WholeBlood=WholeBlood.y)

write.table(comboTW,file= my.dir %&% "GTEx_Tissue-Wide_local_h2_se.txt",quote=F,row.names=F,sep="\t")
write.table(comboTS,file=my.dir %&% "GTEx_Tissue-Specific_local_h2_se.txt",quote=F,row.names=F,sep="\t")